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SUMMARY 


PROBLEM 

Provide  the  underwater  optical  systems  designer  with  a  physical  description  of  multiple 
scattering  of  light.  Derive  simple  approximate  expressions  for  the  target  plane  illumin¬ 
ation  which  are  valid  in  the  region  of  present  hardware  viewing  capabilities. 

RESULTS 

An  iterative  method  is  presented  which  yields  an  exponential  model  for  underwater 
multiple  scattering  in  the  forward  direction.  An  easily  evaluated  effective  attenuation 
coefficient  is  derived  which  al.ows  rapid  prediction  of  flux  through  an  aperture  and 
on-axis  irradiancc.  Extensive  comparison  with  experimental  and  Monte  Carlo  results 
indicate  that  our  simple  approximate  expressions  retain  predictive  value  out  to  seven 
to  ten  attenuation  lengths. 

RECOMMENDATIONS 

Further  studies  into  delimiting  the  range  of  validity,  relaxing  some  approximations  and 
providing  more  rigorous  justification  for  others,  could  improve  the  predictive  utility  of 
the  present  model. 
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INTRODUCTION 


Man's  increased  presence  in  the  underwater  environment  has  led  to  a  need  to  increase 
the  range  of  his  sensory  perceptions  in  this  unfamiliar  milieu.  The  Naval  Undersea 
Center  has  been  particularly  active  in  this  area,  realizing  that  the  ability  ‘o  perform 
useful  tasks  underwater  is  determined  to  a  great  extent  by  the  quantity  and  quality  of 
sensory  inputs.  Recent  years  have  seen  great  improvements  in  our  ability  to  utilize  and 
exploit  the  information  contained  in  optical  and  acoustic  signals. 

In  the  field  of  underwatei  optics,  many  techniques  have  been  introduced  to  enable 
viewing  at  greater  distances.  Though  each  has  its  disadvantages,  the  techniques  of 
range-gating,  volume  scanning,  and  adaptive  scanning  have  yielded  imagery  at  greater 
distances  than  heretofore  possible.  Measurements  of  water  parameters  have  led  to  an 
increased  appreciation  of  the  viewing  problem  and  indicate  that  further  improvements  in 
viewing  range  are  possible. 

Increased  optical  hardware  capability  presents  a  mixed  blessing  to  the  systems  designer. 
To  advantage  is  that  he  can  now  specify  systems  which  are  certain  to  yield  increased 
viewing  ranges  over  prior  techniques.  However,  he  has  difficulty  predicting  the  per¬ 
formance  of  these  systems  since  they  operate  deep  in  the  multiple-scattering  domain. 
Because  the  light  is  scattered  many  times  in  its  travels  from  source  to  target  to  receiver, 
irradiance  values  become  difficult  to  predict.  The  use  ot  conventional  single-scattering 
formulas  can  lead  to  predictions  which  are  in  error  by  orders  of  magnitude. 

The  systems  designer  is  thus  faced  with  the  problem  of  designing  an  underwater  view¬ 
ing  system  which  will  perform  in  the  complicated  multiple-scattering  domain.  Often 
compounding  the  problem  is  the  designer’s  lack  of  familiarity  with  multiple  scattering 
and  the  mathematical  tools  necessary  to  interpret  the  vast  literature  on  this  subject. 

Lack  of  time  and  money  resources,  especially  in  the  preliminary  design  stages,  might 
also  prevent  implementation  by  the  designer  of  Monte  Carlo  or  other  computer-based 
numerical  methods. 

Recognizing  the  system  designer’s  dilemma,  the  Naval  Undersea  Center  (NUC)  has 
expended  considerable  effort  to  give  him  some  methods  that  are  applicable  to  problems 
in  extended  range  viewing.  A  handbook  lias  been  published  (reference  1)  which  uses 
a  hybrid  Monte  Carlo  model  and  the  concept  of  effective  attenuation  coefficients  to 
obtain  a  description  of  image  quality.  Also,  a  canonical  design  procedure,  interrelating 
the  various  system  and  component  parameters,  is  presented.  Another  NUC  report 
(reference  2)  documents  the  Monte  Carlo  program  used  in  the  handbook  and,  by 
listing  the  program,  allows  the  user  to  generate  effective  attenuation  coefficients  for 
cross  sections  not  considered  in  the  handbook. 

The  work  herein  endeavors  to  aid  the  systems  designer  in  two  ways  not  included  in  the 
above  reports.  First,  Part  1  of  this  work  presents  a  simple  analytic  model  of  multiple 
scattering  of  light.  It  is  hoped  that  this  model  will  increase  the  designer’s  appreciation 
of  the  effects  of  the  medium  on  the  system’s  viewing  capability.  Second,  Part  I  of  this 
report  provides  simple  exponential  expressions  for  the  flux  through  an  aperture  and  for 
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the  on-axis  irradianee.  These  expressions  are  so  easy  to  evaluate  that  *’ie  implications 
of  different  scattering  cross  sections,  scattering  coefficients,  and  scattering  to  absorption 
coefficient  ratios  can  in  turn  be  quickly  and  easily  evaluated.  A  further  advantage 
is  that  the  functional  dependence  of  the  effective  attenuation  coefficient  on  the  medium’s 
inherent  optical  properties  is  explicitly  displayed. 

Part  11  of  this  report  exhaustively  compares  the  results  of  Part  I  with  Monte  Carlo 
simulations  and  c  'erimental  data.  Part  II  shows,  for  the  chosen  scattering  cross 
section,  that  all  ti  lee  methods  give  results  in  good  agreement  for  ranges  up  to  six 
scattering  lengths  (i.e.,  generally  seven  to  ten  attenuation  lengths^,  but  that  the  simple 
exponential  model  of  Part  I  breaks  down  at  longer  ranges.  Since  almost  all  state- 
of-the-art  viewing  systems  are  limited  to  six  scattering  lengths  or  less,  the  failure  of 
the  model  at  longer  ranges  will  probaliy  be  of  little  concern  to  all  but  the  mos": 
advanced  and  highest  powered  viewing  systems  applications. 
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INTRODUCTION 


The  p'^blein  of  evaluating  the  intensity  of  a  beam  which  has  passed  through  a  multiple 
scattering  environment  has  occurred  repeatedly  in  widely  varying  areas.  The  passage  of 
charged  particles  through  thick  foils,  the  shielding  and  moderating  of  neutrons,  the 
scattering  of  high  frequency  sound  by  ship  wakes,  and  the  scattering  of  light  both  in 
the  atmosphere  and  in  water  are  all  governed  by  multiple-scattering  effects  if  the 
distances  of  interest  are  great  enough.  Because  of  the  importance  of  multiple  scat¬ 
tering  in  many  fields  and  over  a  long  time  span,  the  scientific  literature  is  replete  with 
contributions  on  the  subject.  In  contrast  with  most  previous  works,  this  report 
endeavors,  at  the  sacrifice  of  some  mathematical  sophistication,  to  provide  a  rather 
intuitive  approach  and  results  which  are  computationally  convenient. 

The  theoretical  approach  which  will  be  used  here  is  known  as  an  iterative  procedure. 
Preisendorfer  (reference  3)  has  called  such  a  procedure  a  natural  tool  .  .  because 
the  light  field  N  in  any  medium  may  be  thought  of  as  the  superposition  of  com¬ 
ponent  fields  N(,)  consisting  of  radiant  energy  having  undergone  i  scatterings,  i  =  1, 

2,  3,  .  .  .,  after  entering  the  medium”.  In  addition,  this  approach  has  other  advan¬ 
tages  as  far  as  the  underwater  systems  designer  is  concerned. 

The  iterative  method  as  used  here  is  strictly  an  integral  approach,  that  is,  it  sums  up 
the  intensities  due  to  rays  traveling  from  the  source  plane  to  the  target  plane.  Since 
there  is  no  need  to  solve  differential  equations,  no  knowledge  of  this  subject  is 
required.  All  the  quantifies  needed  in  the  calculation  are  directly  physically  measurable; 
no  artificial  constructs  such  as  scattering  kernels  or  transmittance  operators  are  intro¬ 
duced.  Also,  the  iterative  equations  are  subject  to  a  simple  implementation  of  the 
small  angle  and  of  the  no-return  backscatter  approximation,  which  are  a  consequence 
of  the  highly  peaked  nature  of  the  volume  scattering  function.  Without  any  evaluation 
of  integrals,  the  resulting  expression  for  the  irradiance  yields  a  straightforward  physical 
interpretafi  m.  Finally,  the  iterative  procedure  used  here  should  have  ap-  .1  to  those 
readers  who  may  be  interested  in  Monte  Carlo  techniques,  since  there  a  c  many 
similarities  in  the  two  methods  of  obtaining  intensities  in  terms  of  the  ray  patiis. 

In  presenting  the  derivation,  the  author  has  tried  to  catalog  ail  approximations  as 
explicitly  as  possible.  Physical  and/or  mathematical  justifications  are  presented  with 
the  hope  that  the  reader  will  be  able  to  assess  the  applicability  of  the  present  results 
to  his  own  particular  problem. 


ANALYSIS 


STATEMENT  OF  PROBLEM 

The  geometry  of  the  problem  is  shown  in  figure  1.  Here,  there  is  a  source  plane 
which  contains  an  initially  specified  distribution  of  light.  At  a  distance  R  from  the 
source  plane  is  the  target  plane  on  which  we  wish  to  know  the  light  distribution.  The 
space  between  these  planes  is  filled  with  water  of  known  optical  properties.  The 
problem  is  to  find  the  distribution  of  the  light  on  the  target  plane  as  a  function  of 
the  distribution  on  the  source  plane,  the  optical  properties  of  the  water,  and  the 
distance  R. 


Two  important  restrictions  allow  us  to  simplify  our  problem  immediately.  First,  we 
shall  neglect  diffractive  effects.  This  is  a  reasonable  approximation  for  most  under¬ 
water  applications  and  is  discussed  in  reference  4.  Second,  we  shall  deal  only  with 
light  intensities  and  not  amplitudes,  i.e.,  we  assume  an  incoherent  process. 

In  view  of  these  restrictions,  we  can  view  light  propagation  in  a  very  naive  way.  The 
source  plane  is  a  specified  source  of  photons,*  each  of  which  propagates  in  a  straight 
line  through  the  medium  until  1)  it  is  absorbed  and  stops,  2)  it  is  scattered  and 
(possibly)  changes  direction  and  continues,  or  3)  it  reaches  the  target  plane  and  stops. 
Figure  1  shows  a  typical  photon  which  undergoes  two  scatterings  before  reaching  the 
target  plane. 

We  define  the  light  distribution  and  water  parameters  in  terms  of  photons  although 
they  could  be  defined  as  easily  in  terms  of  energy  flow.  One  can  hop  back  and 
forth  from  one  description  to  another  through  the  use  of  Planck’s  expression  for  the 
energy  of  a  photon. 


The  distribution  of  light  on  the  source  plane  will  be  described  in  terms  of  the  source 
plane’s  radiance,  NQ  (xQ,  yQ,  0,  ,  <*>,),  defined  so  that  NQ  (x0,  yQ,  0,,  ^,)dx0dy0d  ft, 
is  the  number  of  photons  per  second  leaving  the  area  element  dx()  dy0  centered  at 
(x0,  y0)  and  beading  into  the  solid  angle  dft,,  centered  at  the  angles  0,  and  ^ .  The 
symbols  xQ,  yQ  refer  to  Cartesian  coordinates  and  0  ,  and  sP,  to  polar  coordinates. 

The  total  flux,  FQ,  in  photons  per  second  leaving  the  source  plane  will  be 


(xo’ V°i-^i>dxodni’ 


where  the  limits  of  integration  arc  over  all  xj  and  all  ft,. 


(1) 


*  Photons  in  the  sense  of  Newton's  corpuscles,  i.e ,  photons  which  contain  energy  hut  no  phase  information. 
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Figure  1.  Scattering  Geometry. 


The  light  distribution  at  the  target  is  described  by  the  irradiance  H(xT,  yT  ),  where 
H(xt,  yT)dxT  dyT  is  the  number  of  photons  per  second  arriving  at  the  infinitesimal 
area  dx  T  dyT  centered  at  (x^,  yT).  The  flux,  F,  in  photons  per  second  entering  a 
finite  area  (or  aperture)  A  of  the  target  plane  is 


F=  J  H(xT,yT)dxf.  (2) 

A 

The  basic  optical  properties  that  we  need  are  the  water’s  absorption  coefficient  and  its 
volume  scattering  function.  The  absorption  coefficient,  a,  is  defined  such  that  a*dC 
is  the  probability  of  a  photon  being  absorbed  after  traversing  an  infinitesimal  distance 
d2.  The  volume  scattering  function,  o(6),  is  defined  schematically  in  figure  2  and 
algebraically  by 


o(6)  = 


Ml 

HdV  ’ 


(3) 


where  J(0)df2  is  the  number  of  photons  per  second  scattered  into  a  solid  angle  dfi 
centered  at  0,  and  H  is  the  number  of  photons  per  second  per  unit  area  (i.e.,  the 
irradiance)  entering  the  scattering  volume  dV.  Integrating  equation  3  over  S2  and 
using  equation  2,  we  obtain 


'  IMA 


0 


o(0)sin0d0 


dC, 


since 


the  number  of  photons  per 
the  number  of  photons  per 
dA  •  dC 


second  incoming  =  HdA 
second  scattered  =  J  J(0)di7 


(4) 


(5) 


Because  F  t/Fta  is  just  the  probability  of  being  scattered  in  the  distance  d£,  we 
define  the  scattering  coefficient  s  such  that  the  probability  of  being  scattered  in  the 
distance  af  is 


/• 

/* 

sdC  = 

2nJ 

f  o(0)sin0d0 

dC  = 

j  o(0)dn 

We  can  also  define  the  total  attenuation  coefficient,  a,  as 

a  =  a  +  s, 


(6) 


(7) 


J -..‘Vilts  da  A Mat  -Eft-*  -rvri  >  -  >  ai-V. 


8 


so  that  a  •  dC  is  the  probability  of  either  being  scattered  or  absorbed  in  the  distance 
dC.  Finally,  we  can  define  the  backscatter  coefficient  b  such  that  bsd£  is  the  prob¬ 
ability  of  scattering  in  the  distance  dt  by  an  angle  0  where  ni 2<  0  <  n: 

n 

b  =  2*-  J  o(0)sinOdO.  (8) 

ff/2 

Having  defined  the  problem  and  the  quantities  of  interest,  we  shall  now  outline  the 
procedure  to  be  used  in  arriving  at  a  solution. 

OUTLINE  OF  THE  ITERATIVE  METHOD 

The  fundamental  concept  of  the  iterative  method  is  that  the  irradiance  on  the  object 
plane,  H,  can  be  regarded  as  being  formed  from  the  infinite  series  in  equation  9  where 
H:  is  the  irradiance  contribution  from  all  photons  which  have  undergone  exactly  i 
scatterings. 

00 

H  =  2  H:.  (9) 

i=0 

This  decomposition  is  illustrated  in  figure  3  for  i  =  0,  1,2.  Here,  the  source  plane 
is  labeled  by  O-coordinates,  the  target  plane  for  the  case  of  N  scatterings  is  labeled 
by  x*N+| ,  and  a  set  of  N  new  planes,  called  scattering  planes,  is  introduced  and 
labeled  by  the  coordinates  i  =  1,2,  •••N,  where  N  is  the  number  of  scatterings  for 
the  geometry  when  the  photons  scatter  exactly  N  times.  All  the  planes  are  parallel  and 
the  origins  of  the  planes’  Cartesian  coordinate  systems  Jx;,  yj|  ;)re  all  collinear  with 
a  perpendicular  joining  all  the  planes.  The  distance  between  plane  i  and  i  +  1  is 
denoted  by  rj+1  .  Also,  eacu  plane  (except  the  target  plane)  has  a  set  of  polar  coordin¬ 
ates  assigned  tha:  will  be  used  to  reference  the  direction  from  which  the  photons 
leave  the  plane.  The  polar  coordinates  |  will  be  assigned  to  the  i  -  1  plane. 

The  geometry  just  introduced  allows  us  to  describe  the  iterative  method  in  a  succinct 
fashion.  First,  equation  9  enables  us  to  consider  just  those  photons  which  have 
undergone  exactly  N  scatterings.  These  photons  will  be  broken  down  into  groups, 
called  configurations.  Each  configuration  consists  of  all  photons  which  1)  scatter 
exactly  N  times  and  2)  scatter  within  a  distance  dr  of  r.  =  a[  and  r2  =  a2  and 
•••  rN  =  aN  and  nowhere  else.  The  no-retum  backscatter  approximation  discussed  in 
the  next  section  will  allow  us  to  write  a  simple  expression  for  the  target  plane 
irradiance  of  each  configuration.  Summing  (actually  integrating)  over  the  rj ’s  then 
account  for  all  photons  which  have  scattered  exactly  N  times,  i.c.,  summing  over  all  the 
configurations  yields  the  irradiance  due  to  exactly  N  scatterings.  HN .  Finally,  the 
total  irradiance  on  the  target  plane  due  to  all  the  photons  leaving  the  source  plane 
will  be  obtained  by  summing  the  HN’s,  N  =  0.  I,  •••oc  ,  according  to  equation  9. 

The  iterative  method  thus  involves  segregating  the  light  into  small  components,  solving 
the  problem  for  each  component,  and  then  obtaining  the  total  solution  by  summing 
the  component  solutions.  We  can  do  this  because  of  the  linearity  of  the  processes 
governing  the  transfer  of  radiance. 
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DERIVATION  OF  THE  OBJECT  PLANE  LIGHT  DISTRIBUTION 


THE  NO-RETURN  BACKSCATTER  APPROXIMATION  -  Figure  3(c)  illustrates  two 
fundamentally  different  types  of  rays  (i.e.,  photon  paths)  which  can  occur  in  a  given 
configuration.  Ray  one  has  the  property  of  <  w/2  for  all  i  in  the  configuration, 
i.e.  a  ray  such  as  ray  one  will  never  acquire  a  component  of  motion  directed  from 
the  target  piane  to  the  source  plane.  As  a  consequence,  all  rays  having  this  property 
also  have  the  following  properties:  1)  all  the  photons  enter  each  plane  at  which  they 
are  scattered  from  below  (i.e.,  from  a  direction  towards  the  source  plane)  and  2)  each 
scattering  occurs  sequentially,  i.e..  the  ray  scatters  from  the  first  scattering  plane,  then 
the  second  and  so  on.  Rays  such  as  ray  two  have  the  property  of  0 j  >  ff/2  for 
some  i.  For  a  given  configuration,  this  implies  that  1)  type  two  rays  can  enter  any 
scattering  plane  from  either  above  or  below  and  2)  type  two  rays  can  scatter  from 
the  scattering  planes  in  any  order  (e.g.,  in  figure  3(c)  ray  two  scatters  from  the 
planes  in  the  order  0-2- 1-3).  Clearly,  type  two  rays  undergo  a  more  complex  scattering 
behavior  and  make  our  problem  of  evaluating  the  object  plane  irradiance  more  difficult. 

Fortunately,  in  typical  oceanic  waters  type  two  rays  are  uncommon  compared  to  type 
one  rays.  This  is  a  direct  consequence  of  the  small  value  of  the  backscatter  coef¬ 
ficient  of  ocean  waters.  An  average  of  three  readily  available  backscatter  coefficients 
(references  5  through  7)  indicates  a  value  of  b  =  0.023.  Using  the  results  in  appendix 
A.  we  can  express  the  range  R  at  which  ignoring  type  two  rays  introduces  an  error  e 
in  the  flux  reaching  the  target  plane  as 


R  =  (10) 

b2  c 

where  c  is  the  scattering  to  absorption  coefficient  (s/a)  ratio  and  R  is  given  in  scat¬ 
tering  lengths.  Taking  c  =  5.0  (a  rather  large  value),  e  =  0.05,  and  b  =  0.03,  we 
have  R  =  22  scattering  lengths,  a  distance  much  greater  than  distances  of  interest  here. 
Thus,  we  will  ignore  the  contribution  of  type  two  rays  in  all  the  following  calculations. 

The  remaining  type  one  rays  have  the  important  property  of  interacting  with  each 
scattering  plane  sequentially.  This  means  that  the  radiance  distribution  emerging  from 
the  i  +  1  scattering  plane  is  completely  determined  by  the  radiance  distribution  on 
the  i  plane.  Thus,  to  solve  the  problem  of  the  irradiance  on  the  target  plane  for  a 
configuration  containing  N  scattering  planes,  we  need  only  compute  1)  the  radiance 
distribution  emerging  from  a  scattering  plane  due  to  the  radiance  distribution  on  the 
immediately  preceding  scattering  (or  source)  plane  and  2)  the  irradiance  distribution  on 
the  target  plane  due  to  the  radiance  distribution  on  the  last  scattering  plane. 

TRANSFER  OF  RADIANCE  AND  IRRADIANCE  -  Figure  4  shows  the  geometry  for 
which  we  will  calculate  the  transfer  of  radiance  between  two  planes.  The  geometry 
and  symbols  are  similar  to  those  shown  in  figure  3.  The  known  quantities  are  NQ 
(xQ,  y  ,  0  j,  i ),  the  radiance  distribution  on  the  0-plane,  the  water  parameters,  and 
the  distance  Tj .  It  is  desired  to  obtain  the  outgoing  radiance  of  the  1-plane  in  terms 
of  these  parameters. 
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The  radiant  intensity  emitted  by  a  small  area  dA0  =  dxn  dyQ  of  the  0-plane  centered 
at  (xG,  yQ)  is 


N(xo-  Vo-  *,>dxo  dV 


(11) 


where 


and 


2  i'/z 


0.  =  arc  tan  - 


((x,  -  xQ)2  +  (y ,  -  y0)2] 


(12) 


(13) 


The  total  j  hoton  flux  within  the  solid  angle  dO),  bounded  by  0,  to  0,  +  d0,  and 
to  <p,  +  dp, ,  is  then 


N  ( x  v  0 
'  o'  o'  }o'  : 


p,)dx0  dyG  d 


(14) 


In  free  space  equation  14  would  represent  the  flux  actually  reaching  the  1 -plane.  How¬ 
ever.  because  of  the  scattering  and  absorptive  properties  of  the  water  lying  between  the 
planes,  not  all  of  this  flux  reaches  the  1-plane. 

Considering  the  scattering  interation  first,  the  iterative  model  requires  that  in  any 
given  configuration  we  deal  only  with  the  rays  that  scatter  within  a  distance  dr  of  the 
scattering  plane  (in  this  case  the  1 -plane).  By  considering  the  probability  of  not 
scattering  in  the  distance  r, ,  we  will  find  the  fraction  of  the  flux  (given  by  equation 
14)  which  reaches  the  1 -plane  unscattered. 

If  we  make  the  usual  assumption  that  the  probability  of  scattering  in  a  short  distance 
AC  is  given  by 


sAC,  (15) 

then  the  probability  of  not  scattering  in  that  distance  is 

1  -  sA£  (16) 

and  the  probability  of  not  scattering  in  the  distance  2A£  is 

(1  -  sAC)(  1  -  sA£)  =  (1  -  sA£)2.  (17) 

Now  an  arbitrary  distance,  £, ,  may  be  represented  as 

£,  =  nA£. 
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(18) 


Then,  by  generalizing  equation  17,  we  obtain  the  probability  of  not  scattering 
in  the  distance  C,,  P  (C,  ).  as 

l  ns  i  ■ 


Pns  (C,)  = 


for  n-*-  x  AC-*dC  and 


PnA)=Mm  I 

n-+  x 


(19) 


(20) 


Early  in  this  section  we  stated  that  the  effects  of  absorption  can  be  accounted  for  by 
terminating  absorbed  rays  at  the  point  where  they  are  absorbed.  Although  this  is 
proper,  it  is  inconvenient  to  terminate  rays  in  the  iterative  approach.  Instead,  we 

shall  use  an  equivalent  procedure  (reference  8)  which  allows  each  ray  to  continue  on 

through  but  weights  its  intensity  by  the  probability  of  it  not  being  absorbed  in  the 
distance  traveled.  Since  the  probability  of  being  absorbed  in  the  short  distance  AC 
is  given  as  a  -AC,  we  obtain  by  a  procedure  identical  to  equation  16  through  20.  for 

the  probability  of  not  being  absorbed  in  the  finite  distance  C  ,  P  ((!  ): 

Pn,(V  =  e’aC'-  (21) 


The  probability  of  ocing  neither  scattered  nor  absorbed  in  the  distance  C  is  then 


(  a+s  )C 


-(a+i  ) 


a  r, 


(22) 


The  amount  of  flux  leaving  the  1-plane  from  the  element  dxQ  dy(>  in  the  direction  dfy 
which  arrives  at  the  2-plane,  is  then 

«r  j 

cos  0 

e  N0(xf).  yQ,  01,^i)dxodyodS2].  (23) 


Now  from  figure  4, 


dS2[ 


d A ,  _  dA’  cos2  0, 


(24) 
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Using  equation  23  and  24  we  obtain  for  the  irradiance  on  the  frontal  plane  dA',  of 
the  scattering  volume  dV 


a  r. 


cos  d  cos2 

H(dA,)  =  e  N0(x0,  yc,  0,,  ^,)dx0  dyQ  — ; 


(25) 


Inserting  this  value  into  equation  3,  we  can  calculate  dJ,(x,,  y,,  02,  <p2),  the  scat¬ 
tered  flux  leaving  the  volume  element  dV  on  the  1 -plane  due  to  the  flux  arriving  at 
dV  from  the  vicinity  of  xQ,  yQ,  as 


Otr 


dJ, (x, ,  y,,  02,  <p2)  =  e 


cos  0 ,  COS"  0  . 

1  X!  f.,  ..  a  \  ^  ,U,  ___ - L o  (02i)dV, 


No(xo>  Vo’0!’^^  dxo  dVo 


1 


(26) 


where  02  ,  is  the  angle  between  the  incoming  and  scattered  photon,  i.e.. 


0,,  =  arc  cos[sin0,  cost px  sin02  costp2 

+  sin0,  simp,  sin02  sinv?2  +  cos0,  cos02).  (27) 

At  this  point  it  is  convenient  to  introduce  the  scattering  probability  density  p(0)  such  that 


p(0)  =  -2-M. 


so  that 


(28) 


/■ 


p(0)d£2  = 


(29) 


Noting  that  dV  =  dA,  dr, ,  we  can  now  write  (using  equations  26  and  28) 


Or, 


sdr  cos  6. 

dJ,(x,,  y,,  02,tp2)  =  dA,  — —  e  N0(x0,  yG,  0, ,  *,) 

ri 

cos2  0,  p(02 , )  dx0  dyQ. 


(30) 


The  radiance,  dN,(x,,  y, ,  02,  <p2),  emerging  from  (x, ,  y,)  and  heading  in  the 
direction  (02,  <p2l  due  to  photons  entering  the  scattering  volume  dV  from  an  elemental 
area  surrounding  (xQ,  yQ)  is 


16 


(31) 


dN  j  (x , ,  y , ,  0  2 ,  <p2)  =  c  c°s0*  N0(x0,  yQ,  0,,  v>,) 

ri 

cos20,  A02  ,  )dxQ(  dyQ. 

Integrating  equation  31  over  the  0-plane  coordinates  we  obtain  the  radiance,  N,(x,, 
y,,  0  ,  (£2),  emerging  from  the  point  (x,,  y,)  and  heading  in  the  direction  (02 , 


N,(x,,  y,,  62,  *2)  -  -3 


?// 


cos  0 

e  ^o^xo’  yQ> 


Noting  that 


cos2  0,p(021)dxo  dyQ. 


dn  =  cos3  0idxo  dyQ 


we  arrive  at 


N,  (x i ,  y, ,  02,  *2)  =  sdr 


i 

f  cos  0 

/  e. _ N 

1  J  cos 0,  o 


p(02 ,  )dn, , 


which  achieves  our  desired  result  of  expressing  the  1 -plane  radiance  in  terms  of  tlv* 
0-plane  radiance.  Since  our  derivation  did  not  depend  on  the  detailed  properties  o. 
the  0-plane  or  1 -plane,  the  subscripts  0  and  1  are  dummy  indices  and  we  can  immed¬ 
iately  write 


N2(x2,  y2,  0,,  *>,)  =  sdr. 


^  NlP(032)dS22 


Substituting  the  results  of  equation  34  into  equation  35  we  obtain 


Nj(x2,  y2,  03,  *3)  =  s2dr2  dr, 


COS0j  COS0,~ 


NoP(02.)P(03  2)dnidn2 


Clearly,  by  continuing  this  procedure,  we  can  obtain  the  radiance  emitted  by  the 
NIh  scattering  plane  as 


N 

VXN’  >V  0N+1'  *n+i>  =  n  sdrj 

i+  I 

N 

where  FI  a  indicates  the  product  of  all  the  a 's,  i  =  I,  2.  •••  N,  and  where  the 
!+  1  ‘  ‘ 

quantity  in  brackets  acts  as  an  operator  on  NQ.  Equation  37  thus  provides  us  with 
an  expression  for  the  radiance  leaving  the  last  scattering  plane  of  an  N-seattering  plane 
configuration  in  terms  of  the  initial  radiance,  NL.  and  the  water  characteristics. 

To  complete  the  solution  for  a  given  configuration  we  need  the  target  plane  irradiance 
in  terms  of  the  radiance  leaving  the  N,h  scattering  plane.  The  geometry  for  this 
process  is  shown  in  figure  5.  The  flux.  d2F(xN.  yN,,0N  +  .,  £N  +  1).  leaving  an  area 
surrounding  the  point  (xN,.  yN )  and  heading  into  the  element  of  solid  angle  centered 
at  (0N  +  ] ,  <pN  +  ,  )  is  expressed  according  to  equation  23: 


°rN  + 1 
cos  0 


d  HxN  .  yN .  0N+ t .  v?N  + !  >  e  N+l  Nn (xn  ,  yN .  <?N  + , .  <pN  f , ) 

dxN  dj' v  dr2N  +  j , 


(38) 


where  Nn(xn<  yN ,  0N  +  ,.  ^N+1)  is  the  radiance  leaving  the  N’h  scattering  plane. 
Since 


dS2 


N  +  I 


HA'  cos  ^0 

N  +  1  N  +  1 


*N+  1 


(39) 


and 


(i A  N  +|  +  i  1*xn+  i  dyN+ 1  • 


(40) 


we  have 


a  r 


N  +  1 


l2l,  a  ,  C0S^N  +  i  xi  ,iA  /  d_x, ,dyN_cos%  +  ) 

d  Kxn  ,  yN  .  0N  + 1  '  ^N+ I  ^  “  C  1  n“n*I  (  *  -2 

1  N  +  1 


(41) 


Letting  •7f‘N  be  the  irradiance  on  the  target  plane  due  to  a  configuration  containing  N 
scattering  planes,  we  have 
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(42) 


dJf...  _  ,  C0S^N+t  M  (v  V  0  0 

N  -  t  r''NVXN'  yN’  WN  +  1  ’  ^N+l 


and  the  irradiance  of  the  given  configuration  is 


ttrN+  l 


*N  = 


'N  +  1 


Nn(xN:  yN,  0N  +  1.  ^N  +  i^^N+r 


n 


N  +  1 


(43) 


To  obtain  the  total  target  plane  irradiance,  HN ,  of  all  configurations  which  involve 
exactly  N  scatterings,  we  must  sum  over  all  configurations,  i.e.,  we  must  integrate 
5Cn  over  all  the  do,  i  =  1  ••*N.  For  example. 


for  H 1 


for  H2 


for  H3 


and,  in  general. 


r,=R  r2=R-r3  r,=R-r3-r2 


"'■/  I  f 

r  =0  r  =0  r,  =0 


X.. 


(44) 


(45) 


(46) 


(47) 
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with  the  conventions 


N 

I  r  =  0 


N  +  1 


J 


0 

n  a-  =  1 

i+1 


and  where  the  quantity  in  brackets  in  equation  47  operates  on  r7CN. 
Inserting  the  results  of  equations  37  and  43  into  equation  47  yields  HN  as 


(48) 


Hn  = 


/ 

N 

Or.  1 

1 

r 

r:  =R-  I  r. 

N+l  rl  i+1  J  / 

l 

COS0j 

"i  /  J 

ri  =0  R; 

e  f.  p(0i+l.  i>dr*idri 

N° 

- 

1 

(49) 


where  the  factor  within  the  square  brackets  acts  as  an  operator  and  the  following 
converiions  apply:* 

ffa'K  + 1  )drN  + 1  ^  rN  +  i  ) 

P(0N  +  2’  N+  1  ^  1 

C 


For  example, 


and 


I  r.  =0  for  k  >  £ 
i=k 


f.  = 


cosOj  i  <  N 
1  i  =  N+l 


J 


».■  / 


Or 


COS0, 

N0e  1  dtt, 


"■=S/  If 


a  i 


JL_  a  (R-r. ) 


COS 


0. 


(50) 


(51) 


g  1  cos  (L 

e  Nop(02l)<midSl2dTl.  (52) 


cosO 


0  j 


*  where  0(rN+()  is  an  arbitrary  function  of  rN+1 
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We  achieve  our  goal  of  expressing  the  target  plane  irradiance  in  terms  of  the  initial 
radiance  and  the  water  parameters  by  summing  equation  49  according  to  equation  9: 


In  some  cases  we  are  interested  in  the  flux,  F,  through  an  aperture,  A,  on  the  target 
plane.  This  is  a  somewhat  cruder  measure  of  the  target  light  distribution  and  can  be 
expressed  as 

F  =  /  HdicjT,  541 

A 

where  H  is  given  by  equation  53.  Using  that  equation  and  interchanging  summation 
and  integration  we  obtain 


From  equation  55  we  see  that  the  flux,  F,  also  consists  of  components,  each  com¬ 
ponent  resulting  from  photons  which  have  scattered  exactly  N  times. 


SMALL  ANGLE  APPROXIMATION  -  Although  equation  53  provides  a  complete  solu¬ 
tion  for  the  target  plane  irradiance,  it  suffers  from  two  major  difficulties.  Being  an  infinite 
series,  the  N,th  term  of  which  involves  a  2N-tuple  integral,  equation  53 -poses  serious 
computational  problems.  Second,  the  individual  terms  of  the  series  are  not  in  a  form 
which  is  subject  to  the  most  fruitful  physical  interpretation.  When  the  small  angle 
approximation  is  employed,  equation  53  can  be  cast  in  a  form  which  is  very  revealing. 

The  small  angle  approximation  involves  the  assumption  that  cosff  ~  1  for  all  i.  For 
moderate  ranges  and  for  angular  aperture*'  6  (as  viewed  from  the  source  plane)  which 
fulfill  cos0  ~  1,  this  approximation  will  be  valid  for  the  overwhelming  majority  of 
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the  rays  due  to  the  highly  peaked  nature  of  oceanic  scattering  cross  sections.  For  any 
range  and  ray,  if  0  is  large,  6i  will  clearly  be  large  for  some  i,  and  the  approximation 
does  not  describe  the  target  irradiance  at  these  angles.  Also,  as  the  range  between 
target  and  source  planes  increases,  the  maximum  angle  0max  at  which  the  approxi¬ 
mation  holds  becomes  smaller.  Because  optical  viewing  systems  often  perform  in  the 
region  of  moderate  0  and  R.  we  shall  employ  the  small  angle  approximation  in  all 
that  follows. 


Setting  cosfl|  =  1  in  equation  53  and  noting  that  E  r;  =  R,  we  obtain 

i=l 


N 

f  R-  2  r  i  r 

x  N+l  /  i+1  J  / 

H=  2  sNe'aR  n  /  / 

N=0  i=l  J  J 


P^.,.)d«idr  N0. 


From  appendix  B  we  have  the  formula 


7  f 

i=l  J 


N 

R-  2 
i+1 


Multiplying  equation  56  by  the  right-hand  side  of  equation  57  and  dividing  by  the 
left-hand  side,  we  obtain 


N+l 

n 

i=l 
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Equation  58  is  a  simple  and  readily  interpretable  form  for  the  target  plane  irradiance. 
Imagine  the  target  and  source  plane  separated  by  vacuum  rather  than  water.  Into 
this  vacuum  we  insert  N  scattering  planes.  Each  scattering  plane  has  the  property  of 
unity  probability  of  scattering  a  photon  with  an  angular  probability  density  p(0j+1  j) 
The  planes  are  located  at  rt  =  ar  r2  =  a2  •••,  rN  =  aN .  The  fargct  plane  irradiance 
due  to  the  source  plane’s  radiance  propagating  through  these  scattering  planes  is  then 


^hN^>  is  the  value  of  the  target  plane  irradiance  averaged  over  all  possible  values  of 
the  r.  s.  From  equation  60  we  note  that  this  average  takes  place  without  a  weighting 
factor  dependent  on  any  of  the  r^s.  This  means  that,  given  the  fact  that  N  scat¬ 
terings  occur,  any  of  the  configurations  of  the  scattering  planes  is  equally  likely. 

According  to  equation  58,  we  next  compute  the  weighted  sum  of  the  <^hN^  s- 
These  weights  are  just  the  probability  of  scattering  exactly  N  times.  From  equation 
58  we  see  that  this  probability  is  given  by  the  Poisson  distribution  with  a  mean 
number  Of  scatterings  equal  to  sR. 

Finally,  the  target  irradiance  H  is  obtained  from  this  weighted  sum  by  multiplying 
by  e'aR.  All  the  effects  of  absorption  are  accounted  for  by  this  factor  which  repre¬ 
sents  the  fraction  of  flux  that  reaches  the  target  plane  without  being  absorbed. 

Tile  simplicity  of  equation  58  is  due  to  the  small  angle  and  no-return  backscatter 
approximations.  Unfortunately,  except  for  the  first  two  terms,  it  cannot  be  directly 
evaluated.  In  the  next  section  we  shall  establish  some  preliminary  relations  between 
the  flux  through  an  aperture  and  on-axis  irradiance,  leading  to  approximations  for 
these  quantities  which  can  be  easily  evaluated. 

RELATIONSHIP  BETWEEN  FLUX  AND  ON-AXIS  IRRADIANCE  -  One  of  the 

factors  that  makes  equation  58  difficult  to  evaluate  is  the  possible  generality  of  the 
source  plane  irradiance  NQ.  This  difficulty  can  be  alleviated  by  restricting  equation  58 
to  the  case  where  NQ  is  a  unidirectional,  unipotent  point  source  located  at  the  origin 
and  symbolized  by  N“ps. 


Ny ps  =  6(cos 0,  -  ,  (61) 

where  NyP*  is  of  unit  strength  since 

/  NyP*  dJ?,ydx^  =  1.  (62) 

The  restriction  of  NQ  to  Nyps  can  be  made  without  loss  of  generality  in  the  small 
angle  approximation.  If  Hups(xT,  yT,  R)  is  the  target  plane  irradiance  resulting  from 
a  source  plane  radiance  NqPS(xq,  y0,  0,,  ^ ),  then  the  target  plane  irradiance 
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H(\t,  yT,  R)  of  an  arbitrary  source  N0(x0,  yQ,  0,,  <P,)  can  be  expressed  as 

H(xt ,  yT,  R)  =  |  JJ Hups(xT  -  x0  -  R0JCOSV?,,  yT  -  y0  -  R0,sin^r  R)  • 

Nq(x0’  yO'0\'  )d  v1^, }  (63> 

in  the  small  angle  approximation.  Equation  63  is  just  a  convolution  between  Hups 
and  NQ,  with  Hups  serving  as  an  impulse  response  function. 

The  function  Hups  has  much  utility  in  the  field  of  underwater  optics  and  is  known 
as  the  beam  spread  function.  Due  to  the  reciprocity  theorem  for  underwater  optics 
(reference  9),  Hups  also  yields  (aside  from  scale  factors)  the  point  spread  function, 
which  is  the  image  plane  irradiance  due  to  light  radiated  by  an  omnidirectional  point 
source.  The  Fourier  transform  of  the  point  spread  function,  known  as  the  modulation 
transfer  function  (MTF),  is  used  to  describe  contrast  degradation  as  a  function  of 
range  and  spatial  frequency. 

Our  main  use  of  Hups  will  be  in  conjunction  with  equation  63  and  for  the  case 
where  NQ  becomes  the  radiance  distribution,  Nc,  due  to  a  point  source  conical 
beam  of  strength  FQ  and  solid  angle  £2C  expressed  as 

7jr6(0  °r  e  nc  ) 

Nc  =  •  (64) 

0  0,,  </>,  4  S2C  I 

The  on-axis  (i.e.,  xT  =  0,  yT  =  0)  irradiance  due  to  this  source  is,  from  equation  63, 
H(0,  O,  R)  =  {  fj' rtups  <'xo  -R0,cos*r  ^o  -RVin*>,>  R>  * 

6  (x^dXodfij  |  .  (65) 

Letting 

x  =  R0,cos^);  y  =  R01siny?1  (66) 

and  dropping  the  dummy  index  1,  we  have 


(67) 
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where 


F(fic 

Equation  67  has  been  previously  derived  by  Funk  (reference  10)  for  the  case  of  a 
spherical  geometry. 

Equation  67  allows  us  to  evaluate  the  on-axis  irradiance  of  conical  beams.  The 
restriction  to  on-axis  irradiance  is  usually  of  little  practical  consequence  since  it  is 
common  to  evaluate  imaging  systems  on  the  basis  of  on-axis  performance  (reference 
II).  Equation  67  may  be  used  to  delimit  the  on-axis  irradiance  of  a  square  beam 
by  considering  the  irradiance  of  the  inscribed  and  circumscribed  cones. 

The  most  important  property  of  the  on-axis  irradiance  is  that  it  is  proportional  to  the 
integral  of  the  “ups"  irradiance  across  the  beamwidth.  The  significance  of  this  result 
is  that  the  spatial  derivative  of  the  flux  (i.e.,  the  irradiance)  can  be  expressed  in  terms 
of  the  flux  alone.  In  view  of  the  approximations  we  are  about  to  make,  this  is  a 
welcome  property  since  a  function  itself  is  often  less  sensitive  to  certain  approximations 
than  its  derivative. 

FLUX  THROUGH  AN  APERTURE  The  primary  goal  of  this  section  is  to  compute 
F(£2)  (equation  68).  Knowledge  of  this  function  will  allow  evaluation  of  the  on-axis 
irradiance  as  well  as  energy  distribution  contours  in  the  target  plane. 

Using  equation  58  to  express  Hups,  we  have 

Hups  _  ,.-aR  S-  tsK;  e 
—  C  l*  —  “j-rj 

N=0  N! 

where  from  equation  59 


Substituting  equation  69  into  equation  68,  we  obtain 


I  =  /  Hups  (x,y,R)dxdy. 


(68) 


F(12)  =  e‘aR 


y  (sR> 
N=0  N! 


,N  e-sR 


(71) 
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with 


f“p-  = 


(72) 


Although  equation  71  is  simple  in  form,  consideration  of  equations  72  and  70  shows 
that  its  evaluation  is  non-trivial.  To  obtain  an  indication  of  how  to  proceed,  we 
rearrange  equation  7 1 : 


F(fi) 


(73) 


It  is  now  clear  that  approximating 


by  a  quantity  of  the  form 


(74) 


will  make  equation  73  summable  (to  a  exponential)  and  also  will 
tion  of  only  two  quantities  (a  and  P)  rather  than  the  infinite  set 
justifications  for  approximating  ^nPS^  by 


require  the  evalua- 
|  fNurs  }  •  Some 


(75) 


will  now  be  given. 

Inspecting  equation  75  shows  that  <^f“ps^>  's  rePr°duced  exactly  for  N  =  0  or  1 
and  thus  the  approximation  is  correct  for  small  distances,  i.e.,  for  sR  «  1. 


Second,  equation  75  is  correct  when  the  scattering  function  p(0)  is  a  delta  function, 
i.e.,  when 


p(0)  -  6<-CO_S^lJJ 
2n 


(76) 


Since  oceanic  scattering  functions  are  highly  peaked,  agreement  of  equation  75  in  the 
limit  of  equation  76  is  reassuring. 
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Finally,  equation  75  is  subject  to  a  simple  physical  interpretation,  it  states  that  the 
amount  of  the  flux  (  <  fJJ,ps  >  )  reaching  the  aperture  after  N  scatterings  equals  the 
original  source  strength  (<  f“ps  >  )  times  a  factor  (<  f“ps  >/  <  f“ps  > ),  which 
equals  the  probability  of  being  scattered  by  less  than  an  angle  J2,  this  probability 
being  raised  to  the  Nth  power.  Thus,  equation  75  assumes  that  the  result  of  N 
scatterings  is  the  same  as  the  result  of  N  consecutive  single  scatterings.  For  these 
scatterings  only  the  flux  within  an  angle  ft  contributes  to  the  source  term  for  the  next 
scattering,  and  the  flux  distribution  resulting  from  any  given  scattering  can  be  consoli¬ 
dated  into  a  point  source.  This  interpretation  seems  reasonable  for  scattering  at 
moderate  ranges  by  a  volume  scattering  funcrian  which  is  highly  peaked.  At  longer 
ranges,  equation  75  must  break  down  since  some  of  the  flux  originally  scattered  out¬ 
side  £2  will  return  and  the  succeeding  source  distributions  will  depart  significantly  from 
point  sources. 

A  rigorous  mathematical  argument  concerning  the  domain  of  validity  of  equation  75 
is  beyond  the  scope  of  this  work.  Our  main  justification  for  this  approximation  lies 
in  Part  11  of  this  report,  in  which  formulas  based  on  this  approximation  show  good 
agreement  with  both  Monte  Carlo  and  experimental  results.  Utilizing  equations  75 
and  73.  we  obtain 


F(S2)  =  <fu0ps)>  c'“R  +  s  <Y;ps>/</oPS>  R.  (77) 

To  proceed,  we  must  evaluate  /f“ps\  and  \foPS^  •  From  equations  61.  70. 
and  72  we  have 

<r“oP’>  *  ff  N'T  <18, 

=  yy^cosB,  -  1  )5(*>,  (78) 

Now, 


5(x£)  =  5(xo15(y0)  =  6(x,  -  ROj cos^,  )5fy j 


ROj  sini^j  ). 


(79) 


Evaluating  the  delta  function  over  dJ^  we  obtain 

<T>  =  /  S(x,)8(y,)d^, 


(80) 


so 

<T>  -  ' 


(81) 
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<T‘> . 


^fups^  _ 


Integrating  over  dfi 


ffff  P<02i  >  Ny  dn,dfl,dxTdr. 
/dr> 


<rJ”> 


.  fff 


p(02  )6(x,  )6(y,  )dfi2d^dr1 


A- 


Writing 


=  dxidy,cos3^,  dx,dy 


and  integrating  over  x^,  we  obtain 


<f7‘>  * 


rfp(e  > 

//  ~^-“dX2dr, 


/  dfl 


„  +  y?  l'/2 

=  arc  tan  - - ? _ 


dx2  dx„ dy 


r2  r: 


=  07d02dv2  s  c’ft2 


and  dropping  the  subscript  2,  we  obtain 


</fuPs^>  =  __ 


A  i 

0  K 


P(0)dS2  dr. 
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where  dA  is  the  half-angle  subtended  by  a  circular  aperture  (or  by  the  intersection 
of  a  conical  beam)  centered  on  the  target  plane  as  viewed  from  the  first  scattering 
plane,  i.e.. 


0 


A 


arc  tan 


(89) 


where  rA  is  the  radius  of  the  aperture. 


We  note  that  the  quantity  in  square  brackets  in  equation  88  is  just  the  distribution 
function  associated  with  the  normalized  cross  section  p(0).  The  distribution  function 
is  the  form  in  which  underwater  cross  sections  are  often  tabulated;  however,  evalua¬ 
tion  of  “p  /  requires  the  value  of  the  distribution  function  to  be  averaged 
over  r( .  In  line  with  our  attempt  to  present  the  simplest  possible  description  of 
multiple  scattering  which  retains  some  predictive  value,  we  shall  now  approximate 


<fr>  by 


,-ups  j 
1  ! 


(90) 


That  is,  we  shall  drop  the  average  over  the  scattering  plane’s  location  in  favor  of 
evaluating  the  function  f“ps  for  the  case  of  the  scattering  plane  located  exactly 
halfway  between  the  source  plane  and  the  target  plane.  This  yields 


r,  =R/2 


p(0)dS2, 


(91) 


where  0'  is  the  polar  half-angle  subtended  by  the  circular  aperture  for  an  on-axis 
observer  located  midway  between  the  source  and  target  plane. 

Using  equations  91,  81,  and  77,  we  obtain  for  Rfi).  the  flux  through  a  circular 
aperture  due  to  a  unidirectional  point  source  of  unit  strength, 


F(S2) 


j-a+  s 
e 


p(0)dfi 


(92) 


Due  to  the  linearity  of  the  scattering  process,  for  a  source  of  strength  F0  equation  88 
becomes 


I 


(93) 


Using  equations  92  and  67,  we  obtain  for  H(0,0.R),  the  on-axis  target  plane  irradianee 
due  to  a  point  conical  beam  of  solid  angle  £2C  and  strength  FQ, 


(94) 


We  note  that  both  the  flux  through  an  aperture  and  the  on-axis  irradianee  are  governed 
by  an  effective  attenuation  coefficient,  ael-f,  defined  by 


(95) 


In  addition  to  being  easily  calculable,  ae(f  also  has  reasonable  limiting  forms.  For 
0'~* 0,  aeff— «.  as  would  be  expected  from  a  narrow  beam  where  almost  all  scattered 
photons  remain  outside  the  beam.  For  0'-* n,  a,ff-*a,  which  corresponds  to  the 
intuitive  statement  that  if  all  the  scattered  light  remains  in  the  beam,  this  light  is 
then  attenuated  by  the  absorption  coefficient. 


CONCLUSION 


Equations  93  and  94  represent  the  main  formulas  of  computational  interest  to  the 
systems  designer.  With  very  little  computational  investment,  they  provide  useful 
expressions  for  flux  and  on-axis  ir/adiance.  Part  II  of  this  report  shows  that  (for 
the  cross  section  tested)  equation  93  retains  predictive  value  out  to  six  scattering 
lengths  or  seven  to  ten  attenuation  lengths,  depending  on  the  s/a  ratio. 

We  have  used  equation  58  as  a  starting  point  for  deriving  expressions  for  on-axis 
irradiance  and  flux  through  an  aperture.  It  also  serves  as  an  excellent  departure  for 
the  derivation  of  the  water’s  MTF.  By  taking  the  Fourier  transform  of  equation  58 
and  applying  the  small  angle  approximation  consistently,  a  series  expression  for  the 
MTF  can  be  obtained.  By  suitably  applying  and  extending  the  model  other  quantities, 
such  as  off-axis  irradiance.  can  also  be  derived. 

The  iterative  method  has  shown  its  ability  to  yield  results  in  agreement  with 
Monte  Carlo,  experimental,  and  analvtic  approaches  and  should  thus  be  considered 
a  useful  and  important  approach  to  multiple  scattering  problems. 
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APPENDIX  A 

THE  NO-RETURN  BACKSCATTER  APPROXIMATION 


We  will  confine  our  attention  to  scattering  on  a  line  since  for  this  case  the  problem 
is  exactly  solvable.  In  figure  A-l,  FQ  is  a  source  of  photons  located  at  x  =  0,  where 
all  these  photons  initially  travel  to  the  right  on  the  line  [0,  °o  ]  and 


u(x) 

v(x) 

b,a,s 

A 

bsA 

aA 


flux  at  point  x  traveling  to  the  right 
flux  at  point  x  traveling  to  the  left 

baekscatter,  absorption,  and  scattering  coefficient,  respectively 
a  small  increment  in  length  along  the  x-axis 
probability  of  being  backscattered  in  the  interval  A 
probability  of  being  absorbed  in  the  interval  A. 


By  considering  the  various  loss  and  gain  processes,  we  can  write  for  u(x) 


u(x  +  A)  =  u(x)  -  aAu(x)  -  bsAu(x)  +  bsAvix  +  A). 


(A- 1) 


In  this  equation,  u(x)  is  the  original  flux,  aAu(x)  is  the  amount  lost  by  absorption, 
bsAu(x)  is  the  amount  lost  by  baekscatter,  and  bsAv(x  +  A)  is  the  amount  of  flux 
gained  by  the  rightward-moving  beam  due  to  baekscatter  of  the  leftward-moving  beam 
(i.e.,  the  returned  baekscatter). 

Noting  that 

v(x  +  A)  ~  v(x)  +  —  A  ,  (A-2) 

dx 

x=x 

we  have,  to  first  order  in  A, 


u(x  +  A)  =  u(x)  -  (a  +  bs)Au(x)  +  bsAv(x)  .  (A-3) 


Similarly,  for  v, 


v(x)  =  v(x  +  A)  -  (a  +  bs)Av(x)  +  bsAu(x)  .  (A-4) 


Dividing  equations  A-3  and  A-4  by  A,  and  rearranging  and  letting  A~K),  we  obtain 


(A-5) 


—  =  -qu(x)  +  bsv(x) 
dx 
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U(x)  - ► 

-  v(x) 


u(0)=  Fo 


Figure  A-l.  One-Dimensional  Baekscatter  Geometry. 
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and 


(A-6) 


—  =  +qv(x)  -  ’osu(x)  , 

where  q  =  a  +  bs. 

Differentiating  equation  A-5  and  using  equations  A-5  and  A-6,  we  obtain 


^  du  +  bs  dy 
dx1'  dx  dx 


The  question  is.  at  what  range  x  =  r  does  the  flux  predicted  by  the  simplified  result, 
equation  A-17,  begin  to  depart  significantly  from  equation  A-13,  which  includes  the 
effects  of  returned  backscatter.  If  we  define  a  fractional  error  e  «  1  as  a  significant 
departure,  then  we  require 


(A- 18) 
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or 

+(-  yj  -  (bs)2  +  q), 

e  =  1  +  e  . 

Letting 

a  =  s/c 
c  =  s/a  ratio 
sr  =  R 

we  have 


(A- 19) 


(A-20) 


+<-l/c  F  +  2  b  c  +  1/c  +  b  >R 

=  1  +  c  (A-21) 


for  ocean  waters 


so  that 


2bc  «  1  , 


(A-22) 


(1  i-  2bc)H  =  1  +  be  - 

2  ' 

Inserting  equation  A-23  into  A-21, 


(A-23) 


b2cR 


=  1  +  C  , 


but  since  e  «  1  wc  obtain 


(A-24) 


(A-25) 


37 


gtff 


APPENDIX  B 

PROOF  OF  EQUATION  57 


Equation  57  asserts  that 


r  f 

i=l  J 


N 

R-  I  r 
i+1 


dr.  = 


R^_ 

N! 


0 


First,  we  will  verify  this  formula  for  N  =  1.  2.  3. 
For  N  -  1. 


1 


rR-  I  r,  r  R 

i /  1+1  jr'=/ 


0 


For  N  =  2, 


R-  r 


R  /•  R-r, 


n 

i-1 


i+1 


dr  = 


dr,  dr2 


0  0 


/'R 

-  J  (R-r,  )dr2 
0 


(""■=) 


R 


=  r2  R_2  =  R : 


For  N  =  3, 


■7 

i=  1  J 


R-  I  r . 


i+1 


J 


R  f  R=r3  f  R-r3  -r, 


dr;  = 


dri  dr2dr3 


0  0  0 
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J  ftffrNtflmfttiiill  '■iUrtf' i  fill 


■jm 


(R  -  r3  -  r2)dr2dr3 


J  0 


R3  =  R3 
6  3!  ‘ 


Now,  if  we  define 


N 

n 

i=l 


fR- 

J 


N 

v 

i+1 


r 

j 


R- 


N  N  +  1 


dr,  •••  drN  ~  FN(R) 


and  assume 


dN 

Fw  (R)  = 

N  N! 


we  need  only  prove  that 


dN  +  1 

Fn  +  .(R)  =  £ - 

N+1  (N+l)! 
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But 


Letting  R  -  rN  +  ,  =  Rf  , 


then 


but  the  quantity  in  brackets  is  FN(R'), 
so 

0 

But 

drN.|  =t!R 
for 

rN+1  =  0  R'  =  R 

rN  +  1  =  R  R'  =  0 

so  that 

fn+,W  =  -/  -C^NdR' 

R 
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and  evaluating  the  integral  yields 


+ 1 (R) 


RN+' 

(N+1)! 
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INTRODUCTION 


The  purpose  of  this  section  is  to  evaluate  the  predictions  of  radiant  energy  flux 
through  an  aperture  developed  from  Gordon's  exponential  model.  This  will  be  done  by 
comparing  them  with  the  predictions  of  a  Monte  Carlo  simulation  model  and  with 
experimental  data. 

The  radiant  energy  flux  passing  through  an  aperture  has  significance  in  underwater 
optics  for  two  reasons.  First,  it  is  important  in  areas  concerning  the  power  distribution 
of  a  collimated  source  (e  g.,  a  laser)  after  the  beam  has  travelled  a  given  distance  under¬ 
water;  for  example,  target  illumination  or  amount  of  energy  reaching  a  receiver. 

Second,  the  flux  through  an  aperture  is  related  to  the  on-axis  irradiance  of  a  conical 
beam,  which  is  significant  because  the  available  irradiance  is  an  important  design  param¬ 
eter  for  underwater  optical  imaging  systems. 

Because  of  the  significance  of  the  radiant  energy  flux,  rather  lengthy  comparisons  will 
be  made  of  this  quantity  as  a  function  of  aperture  angle,  range,  and  scattering  to 
absorption  coefficient  (s/a)  ratio.  It  is  hoped  that  these  comparisons  will  serve  as  a 
guide  to  the  various  rcgions  of  validity  for  the  results  of  Gordon. 
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COMPARISON  SOURCES 


EXPONENTIAL  MULTIPLE-SCATTERING  MODEL 

The  exponential  multiple-scattering  model  developed  by  Gordon  in  Part  I  is  presently 
formulated  to  predict  the  radiant  energy  flux  produced  by  a  unidirectional  point  source 
for  a  given  combination  of  range,  aperture  angle,  and  s/a  ratio. 

To  illustrate  the  geometry  of  this  model,  let  a  unidirectional  point  source  illuminate  a 
screen  placed  perpendicular  to  the  light  source  axis.  If  FQ  is  the  radiant  energy  emitted 
by  the  unidirectional  point  source  and  F  is  the  energy  flux  intercepted  by  the  screen 
placed  at  a  distance  R  inside  a  cone  of  half-angle  7  (a  circular  aperture  of  radius  r), 
then  the  ratio  of  these  will  yield  the  fraction  of  energy  that  remains  inside  the  circular 
aperture.  Tiiis  arrangement  is  shown  in  figure  1.* 

The  basic  equation  of  the  model  as  given  by  Gordon  is 

JL=c  i +  %fi  ,  (i) 

where 

a  =  total  attenuation  coefficient  (meters'1). 

P(a)  5  (o)  (fl)/s  =  probability  density  function  of  volume  scattering. 

R  =  range  (meters). 

7'  =  half-angle  of  the  cone  measured  from  the  midpoint  of  the  range 

(degrees), 

tE-  =  normalized  radiant  energy  flux  through  the  .arget  aperture. 

‘‘o 

It  has  been  established  that  the  ratio  c  =  s/a  is  easily  measured,  but  the  quantity  s 
itself  is  not.  To  write  equation  1  in  terms  of  the  ratio  c,  the  quantity  crR,  which 
defines  the  range  as  measured  in  attenuation  lengths,  must  be  written  in  terms  of  c. 
Substituting  the  definition  of  the  total  attenuation  coefficient 

ot  =  a  +  s,  (2) 

aR  can  be  written  as 

aR  =  (a  +  s)R 

=  aR(  1  +  c).  (3) 


*Figures  are  placed  ct  the  end  of  Part  II. 
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Rewriting  equation  1,  we  obtain 


_F_  _  J-aR  +  sR  f  ,P(o)<m] 

f0  Jy 

=  e[-aR  +  eaR J  ,P(o)dft) 

„  c  hR  +  rff/yHoXin] 

-  e-  4'  '  '"f'/r'H.  (4, 

Front  equation  4  it  can  be  seen  that  the  radiant  energy  flux  for  a  unidirectional  point 
source  decreases  exponentially  with  distance.  It  can  also  be  seen  that  an  effective 
attenuation  coefficient,  a  j-f,  can  be  defined  such  that 

aeff  =  a[  1  "  r^cfy  P(o)d«]  •  (5) 

Single-scatter  underwater  visibility  models  predict  a  normalized  radiant  energy  flux 
given  by 


which  ■  orresponds  to 


F  -  e-aR 
Fo 


fP(o)dS2  =  0. 


Therefore,  the  effective  attenuation  coefficient  aeff  is  bounded  on  one  side  by  a. 
On  the  other  side,  the  effective  attenuation  coefficient  is  bounded  by  a-  s  since 


and 


P(a)dJ2  >  0 


c  _  s. 

1  +  c  a' 


Hence, 


a  >  aeff  >  a. 


(6) 


MONTE  CARLO  SIMULATION  MODEL 


In  order  to  simulate  increasingly  complex  underwater  imaging  systems,  Funk  has 
developed  a  method  for  describing  the  propagation  of  light  underwater  that  includes 
the  effects  of  multiple-scatterings  (references  1  and  2).  In  this  method  Monte  Carlo 
techniques  are  used  to  calculate  effective  attenuation  coefficients  for  illuminating,  image¬ 
forming,  and  backscattered  light;  these  coefficients  permit  the  most  significant  effects 
of  multiple-scattering  to  be  incorporated  into  ar.  underwater  visibility  model. 
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The  simplest  form  of  the  Monte  Carlo  calculations  for  predicting  the  radiant  energy 
flux  from  a  unidirectional  point  source  is  a  specialized  three-dimensional  random  walk. 
The  scattering  and  absorption  properties  of  the  water  are  found  in  terms  of  probability 
distribution  functions  th»t  characterize  the  direction  and  length  of  each  step  in  the 
random  walk. 

DATA  FROM  DUNTLEYS  EXPERIMENTS 

For  a  number  of  years,  S.  Q.  Duntley  of  Scnpps  Institution  of  Oceanography,  Visibility 
Laboratory,  has  been  investigating  the  propagation  of  light  in  water.  During  this 
investigation,  numerous  experiments  were  performed  in  both  fresh  water  at  Diamond 
Island,  New  Hampshire,  and  simulated  ocean  water  of  various  s/a  ratios  in  a  laboratory 
tank.  A  laser  was  used  as  the  light  source,  and  measurements  of  the  radiant  energy 
flux  were  made  at  varying  ranges  and  collector  angles.  The  experimental  procedure  was 
to  measure  the  energy  at  an  off-axis  collector  angle  7'  at  a  radius  r  from  the  source. 
The  data  were  then  integrated  to  obtain  the  total  power  on  a  spherical  cap,  with  7'  as 
the  half  angle  of  the  cap.  These  experiments  are  described  in  reference  3 
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COMPARISON  OF  EXPONENTIAL  MODEL  WITH 
TWO  OTHER  SOURCES 


Using  equation  4,  it  is  now  possible  to  compare  the  predictions  of  the  exponential 
model  with  those  of  the  Monte  Carlo  multiple-scattering  model  and  the  experimental 
data  accumulated  by  Duntley. 

Theoretical  predictions  were  provided  by  Funk  for  the  chosen  cases  from  the  Monte 
Carlo  simulation  model,  and  equation  4  was  programmed  on  a  desk  calculator  to  give 
predictions  from  the  exponential  model.  Representative  sets  were  chosen  from  Duntley’s 
data  to  serve  as  a  basis  from  which  to  assess  the  accuracy  of  the  two  visibility  models. 

The  /yP(o)dS2  was  interpolated  from  the  table  of  the  experimental  volume  scattering 
distribution  function  as  calculated  by  Duntley  (table  1  *).  The  seven  aperture  cone 
half-angles  chosen  by  Funk  for  the  Monte  Carlo  model  were  used  here.  The  angles 
that  Duntley  defines  in  his  experiments  and  those  used  in  Gordon’s  theoretical  model 
are  related  through  the  expression 


tan  y'  =  2  tan  y. 

The  angles  y  (and  7')  and  the  corresponding  scattering  probabilities  are  shown  in 
table  2. 

Tables  3  through  1 1  show  the  radiant  energy  flux  as  a  function  of  aperture  and  range 
for  three  s/a  ratios.  The  results  of  the  calculations  using  equation  4  are  shown  in 
tables  3,  6,  and  9  for  the  ratios  3.56,  2.28,  and  1.48,  respectively.  In  a  similar 
manner,  tables  4,  7,  and  10  display  the  predictions  of  the  Monte  Carlo  simulation 
model,  and  tables  5,  8,  and  1 1  display  the  experimental  results  as  compiled  by  Duntley. 

Figures  2  through  22  show  the  predictions  of  the  exponential  multiple-scattering  model 
(designated  as  E-model  in  the  graphs  for  brevity)  as  straight  lines  on  semi-logarithmic 
paper.  For  comparison,  the  predictions  of  the  Monte  Carlo  simulation  model  (desig¬ 
nated  as  M-model)  and  the  experimental  data  of  Duntley  (designated  as  D-data)  appear 
as  crosses  and  circles,  respectively,  on  the  appropriate  exponential  model  line. 

RELATIVE  ERROR  AS  A  FUNCTION  OF  RANGE 

Figures  2  through  22  indicate  that  both  the  exponential  and  Monte  Carlo  simulation 
model  predictions  give  very  close  agreement  with  the  set  of  selected  experimental  data. 
In  order  to  see  this  agreement  more  easily,  a  simple  type  of  error  analysis  was  per¬ 
formed  to  give  an  idea  of  the  relative  differences  between  the  results  of  the  two 
models  and  the  experimental  data. 


*  Tables  are  placed  at  the  end  of  Part  II. 
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Because  the  exponential  model  of  Gordon  is  based  upon  the  assumption  of  small  angle 
forward  scattering,  the  validity  of  the  small  angle  approximation  is  reduced  as  the 
number  of  scatterings  increases,  which  tends  to  make  the  model  predictions  less  accu¬ 
rate.  To  show  this,  the  unit  of  measurement  for  the  range  was  changed  from  attenua¬ 
tion  length  to  scattering  length.  Now.  the  exponential  model  should  tend  to  break 
down  at  the  same  range  in  all  eases;  i.e.,  the  model  should  break  down  after  the  same 
number  of  scattering  lengths,  regardless  of  the  s/a  ratio.  The  scattering  length,  sR.  can 
be  found  from  the  attenuation  length  aR  through 

sR  =  .  (8) 

1  +  S. 
s 

where 

sR  =  range  (scattering  lengths), 
aR  =  range  (attenuation  lengths), 
a  =  absorption  coefficient, 
s  =  scattering  coefficient. 

To  find  this  average  error,  an  error  equation  of  the  form 


Error  - 


7  angles 

E 


In  ( theoretical  model  flux 
\  experimental  flux 


was  used  to  find  an  average  error  value  over  all  of  the  seven  aperture  cone  half-angles. 
Table  12  shows  the  results  of  these  calculations.  The  errors  between  the  exponential 
model  predictions  and  the  experimental  values  are  denoted  by  E/D.  while  the  errors 
between  the  Monte  Carlo  simulation  model  predictions  and  the  experimental  values  are 
denoted  by  M/D.  Figures  23  through  25  display  these  results  for  the  three  s/a  ratios 
investigated. 

RELATIVE  ERROR  AS  A  FUNCTION  OF  APERTURE  ANGlE 

Another  way  of  examining  the  differences  between  the  predicted  radiant  energy  flux 
and  the  experimental  data  is  to  calculate  the  error  averaged  over  all  the  ranges  for  the 
various  aperture  cone  half-angles.  In  this  case,  the  formula  used  to  calculate  the 
average  error  is  given  by 


10  ranges 


ln  /  theoretical  model  flux 
l  experimental  flux 


Error  = 
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The  results  of  these  calculations  for  the  exponential  model  (E/D)  and  the  Monte  Carlo 
simulation  model  (M/D)  are  summarized  in  table  13.  Figures  26  through  28  display 
these  results  for  the  three  s/a  ratios  investigated. 


CONCLUSIONS 


By  examining  figures  23  through  25,  it  can  be  seen  that  the  differences  between  the 
predictions  of  the  exponential  model  and  the  multiple-scattering  simulation  model  for 
ranges  up  to  six  scattering  lengths  are  not  very  significant.  In  terms  of  attenuation 
lengths,  this  range  corresponds  to  seven  to  ten  attenuation  lengths,  depending  on  the 
s/a  ratio.  The  exponential  model  does  tend  here  to  slightly  better  predict  this  particu¬ 
lar  set  of  experimental  data,  but  it  is  also  apparent  from  these  graphs  that  when  the 
exponential  model  does  break  down  at  the  longer  ranges,  the  Monte  Carlo  simulation 
method  still  continues  tc  offer  reasonable  predictions  for  calculating  effective  attenua¬ 
tion  coefficients. 

Looking  at  figures  26  through  28,  it  is  apparent  that  the  accuracy  of  the  exponential 
model  is  greater  for  small  angles  (less  than  2°)  and  large  angles  (greater  than  10°),  but 
that  the  Monte  Carlo  simulation  model  appears  better  able  to  predict  the  experimental 
results  for  the  intermediate  angles.  Again,  the  differences  are  undoubtedly  not  very 
significant. 

In  addition,  it  must  be  stressed  that  the  accuracy  of  the  two  models  has  been  assessed 
here  on  the  basis  of  only  one  set  of  experimental  data.  While  there  is  no  reason  to 
believe  that  the  data  are  biased  in  any  manner,  such  a  limited  comparison  may  not 
yield  the  same  results  as  one  supported  on  a  broader  base  of  experiment.  Such  a 
comparison  may  show  slight  differences  in  the  relative  accuracy  of  the  two  models. 

Therefore,  in  the  final  analysis,  it  can  be  stated  that  the  ranges  up  to  seven  to  ten 
attenuation  lengths,  this  simple,  exponential  formula  rivals  the  Monte  Carlo  model  in 
accuracy  of  predicting  radiant  energy  flux  or  on-axis  irradiance  due  to  conical  beams 
and  surpasses  it  in  ease  of  use.  When  supplied  with  a  table  of  volume  scattering 
distribution  functions,  one  can  easily  predict  what  fraction  of  the  energy  produced  by 
a  narrow-beam  light  source  will  reach  a  target  of  known  size  located  at  selected 
distances  from  the  source.  The  limitation  to  ranges  of  less  than  seven  to  ten  attenua¬ 
tion  lengths  is  not  a  disabilitating  one  since  the  most  sophisticated  underwater  viewing 
systems  now  available  can  detect  an  object  at  a  maximum  range  of  only  about  seven 
attenuation  lengths.  However,  for  questions  other  than  fiux  through  an  aperture  or 
on-axis  irradiance  or  for  ranges  greater  than  six  scattering  lengths,  one  must  at  present 
look  for  tools  other  than  the  exponential  model,  such  as  Monte  Carlo  analysis. 
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Figure  6.  Radiant  Fnergy  Flux  as  a  Function  of  Range  for  s/a  =  3.56,  7=  10 


Figure  17.  Radiant  Energy  Flux  as  a  Function  of  Range  for  s/a  =  1.48.  7  =  1.00°. 
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Figure  20.  Radiant  Energy  Flux  as  a  Function  of  Range  for  s/a  =  1.48,  y  =  10.00 


Figure  22.  Radiant  Energy  Flux  as  a  Function  of  Range  for  s/a  =  1.48,  7=  100.00°. 
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Figure  25.  Average  Error  as  a  Function  of  Range  for  s/a  =  1.48. 
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T  ible  1.  Volume  Scattering  Distribution  Functions  From  S.  Q.  Duntley. 
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